home *** CD-ROM | disk | FTP | other *** search
/ The Very Best of Atari Inside / The Very Best of Atari Inside 1.iso / mint / mntlb20 / lib / _mulsf3.cpp < prev    next >
Text File  |  1992-03-30  |  4KB  |  159 lines
  1. | mjr: not needed on the TT
  2.  
  3. #ifndef    __M68881__
  4. # ifdef    sfp004
  5.  
  6. | single precision floating point stuff for Atari-gcc using the SFP004
  7. | developed with gas
  8. |
  9. | single floating point multiplication routine
  10. |
  11. | M. Ritzert (mjr at dfg.dbp.de)
  12. |
  13. | 7. Juli 1989
  14. |
  15. | revision 1.1: adapted to the gcc lib patchlevel 58
  16. | 4.10.90
  17.  
  18. comm =     -6
  19. resp =    -16
  20. zahl =      0
  21.  
  22.     .text
  23.     .even
  24.     .globl    __mulsf3, ___mulsf3
  25.  
  26. __mulsf3:
  27. ___mulsf3:
  28.     lea    0xfffa50,a0
  29.     movew    #0x4400,a0@(comm)    | load first argument to fp0
  30.     cmpiw    #0x8900,a0@(resp)    | check
  31.     movel    a7@(4),a0@
  32.     movew    #0x4427,a0@(comm)
  33.     .long    0x0c688900, 0xfff067f8
  34.     movel    a7@(8),a0@
  35.     movew    #0x6400,a0@(comm)    | result to d0
  36.     .long    0x0c688900, 0xfff067f8
  37.     movel    a0@,d0
  38.     rts
  39.  
  40. # else    sfp004
  41.  
  42. | single floating point multiplication routine
  43. |
  44. | written by Kai-Uwe Bloem (I5110401@dbstu1.bitnet).
  45. | Based on a 80x86 floating point packet from comp.os.minix, written by P.Housel
  46. |
  47. |
  48. | Revision 1.2, kub 01-90 :
  49. | added support for denormalized numbers
  50. |
  51. | Revision 1.1, kub 12-89 :
  52. | Created single float version for 68000. Code could be speed up by having
  53. | the accumulator in the 68000 register set ...
  54. |
  55. | Revision 1.0:
  56. | original 8088 code from P.S.Housel for double floats
  57.  
  58. BIAS4    =    0x7F-1
  59.  
  60.     .text
  61.     .even
  62.     .globl    __mulsf3, ___mulsf3
  63.  
  64. __mulsf3:
  65. ___mulsf3:
  66.     lea    sp@(4),a0
  67.     moveml    d2-d5,sp@-
  68.     moveml    a0@,d4/d5    | d4 = v, d5 = u
  69.     subw    #8,sp        | multiplication accumulator
  70.  
  71.     movel    d5,d0        | d0 = u.exp
  72.     swap    d0
  73.     movew    d0,d2        | d2 = u.sign
  74.     lsrw    #7,d0
  75.     andw    #0xff,d0    | kill sign bit
  76.  
  77.     movel    d4,d1        | d1 = v.exp
  78.     swap    d1
  79.     eorw    d1,d2        | d2 = u.sign ^ v.sign (in bit 31)
  80.     lsrw    #7,d1
  81.     andw    #0xff,d1    | kill sign bit
  82.  
  83.     andl    #0x7fffff,d5    | remove exponent from u.mantissa
  84.     tstw    d0        | check for zero exponent - no leading "1"
  85.     beq    0f
  86.     orl    #0x800000,d5    | restore implied leading "1"
  87.     bra    1f
  88. 0:    addw    #1,d0        | "normalize" exponent
  89. 1:    tstl    d5
  90.     beq    retz        | multiplying zero
  91.  
  92.     andl    #0x7fffff,d4    | remove exponent from v.mantissa
  93.     tstw    d1        | check for zero exponent - no leading "1"
  94.     beq    0f
  95.     orl    #0x800000,d4    | restore implied leading "1"
  96.     bra    1f
  97. 0:    addw    #1,d1        | "normalize" exponent
  98. 1:    tstl    d4
  99.     beq    retz        | multiply by zero
  100.  
  101.     addw    d1,d0        | add exponents,
  102.     subw    #BIAS4+16-8,d0    | remove excess bias, acnt for repositioning
  103.  
  104.     clrl    sp@        | initialize 64-bit product to zero
  105.     clrl    sp@(4)
  106.  
  107. | see Knuth, Seminumerical Algorithms, section 4.3. algorithm M
  108.  
  109.     movew    d4,d3
  110.     mulu    d5,d3        | mulitply with bigit from multiplier
  111.     movel    d3,sp@(4)    | store into result
  112.  
  113.     movel    d4,d3
  114.     swap    d3
  115.     mulu    d5,d3
  116.     addl    d3,sp@(2)    | add to result
  117.  
  118.     swap    d5        | [TOP 8 BITS SHOULD BE ZERO !]
  119.  
  120.     movew    d4,d3
  121.     mulu    d5,d3        | mulitply with bigit from multiplier
  122.     addl    d3,sp@(2)    | store into result (no carry can occur here)
  123.  
  124.     movel    d4,d3
  125.     swap    d3
  126.     mulu    d5,d3
  127.     addl    d3,sp@        | add to result
  128.                 | [TOP 16 BITS SHOULD BE ZERO !]
  129.     moveml    sp@(2),d4-d5    | get the 48 valid mantissa bits
  130.     clrw    d5        | (pad to 64)
  131. 2:
  132.     cmpl    #0x0000ffff,d4    | multiply (shift) until
  133.     bhi    3f        |  1 in upper 16 result bits
  134.     cmpw    #9,d0        | give up for denormalized numbers
  135.     ble    3f
  136.     swap    d4        | (we''re getting here only when multiplying
  137.     swap    d5        |  with a denormalized number; there''s an
  138.     movew    d5,d4        |  eventual loss of 4 bits in the rounding
  139.     clrw    d5        |  byte -- what a pity 8-)
  140.     subw    #16,d0        | decrement exponent
  141.     bra    2b
  142. 3:
  143.     movel    d5,d1        | get rounding bits
  144.     roll    #8,d1
  145.     movel    d1,d3        | see if sticky bit should be set
  146.     andl    #0xffffff00,d3
  147.     beq    4f
  148.     orb    #1,d1        | set "sticky bit" if any low-order set
  149. 4:    addw    #8,sp        | remove accumulator from stack
  150.     jmp    norm_sf        | (result in d4)
  151.  
  152. retz:    clrl    d0        | save zero as result
  153.     addw    #8,sp
  154.     moveml    sp@+,d2-d5
  155.     rts            | no normalizing neccessary
  156.  
  157. # endif    sfp004
  158. #endif    __M68881__
  159.